#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/times.h>
#include "mytime.h"

#define N_DEF 200
#define TOL_DEF 0.0000001
#define MAXIT_DEF 500

/****************/
int *despl;
int *count;
/****************/

double
jacobi_llc (double **a, int n, double *b, double *x, double tol, int maxit)
{
  int i, j, it = 0;
  double sum, resi = 10.0, *new_x;

  new_x = (double *) malloc (n * sizeof (double));

  for (i = 0; i < n; i++)
    x[i] = b[i];
  /* while ( (it < maxit) && (resi > tol) ) { */
  while (it < maxit)
    {
      resi = 0.0;
#pragma omp target device(mpi)
#pragma omp parallel for private (sum, j) reduction (+ : resi)
#pragma llc result (&new_x[i], 1)
      for (i = 0; i < n; i++)
	{
	  sum = 0.0;
	  for (j = 0; j < n; j++)
	    sum += a[i][j] * x[j];
	  resi += fabs (sum - b[i]);
	  sum += -a[i][i] * x[i];
	  new_x[i] = (b[i] - sum) / a[i][i];
	}
      for (i = 0; i < n; i++)
	x[i] = new_x[i];
      it++;
    }

  free (new_x);

  return resi;
}

int
main (int argc, char *argv[])
{
  double tol;
  double **a, *b, *x;
  int i, j, n, nfilas, diag, maxit;
  int nprocs, miid;

  CLOCK_TYPE chrono;
  double llc_time;
  double llc_resi;

  nprocs = LLC_NUMPROCESSORS;
  miid = LLC_NAME;

  n = N_DEF;
  maxit = MAXIT_DEF;
  tol = TOL_DEF;

  switch (argc)
    {
    case 4:
      tol = atof (argv[3]);
    case 3:
      maxit = atoi (argv[2]);
    case 2:
      n = atoi (argv[1]);
    }


  /* Por simplicidad supondremos que n es divisible por el numero de procesos: nprocs */
  nfilas = n / nprocs;

  ReservaMatriz (&a, n, n);
  ReservaVector (&b, n);
  ReservaVector (&x, n);

  if (miid == 0)
    {
      GeneraVector (b, n);
    }
  MPI_Bcast (b, n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  /* matriz aleatoria */
  if (miid == 0)
    {
      GeneraMatriz (a, n, n);
      /* Para garantizar la convergencia del metodo de Jacobi */
      /* Convierte en diagonal dominante la matriz distribuida */
      diag = 0;
      for (i = 0; i < n; i++)
	{
	  for (j = 0; j < n; j++)
	    a[i][diag] += a[i][j];
	  diag++;
	}
    }
  /* envio de la matriz */
  for (i = 0; i < n; i++)
    MPI_Bcast (a[i], n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Barrier (MPI_COMM_WORLD);

  /****************************************************************************
   jacobi llc */
  CLOCK_Start (chrono);
  llc_resi = jacobi_llc (a, n, b, x, tol, maxit);
  CLOCK_End (chrono, llc_time);
  MPI_Barrier (MPI_COMM_WORLD);
/*  CompruebaSolucion (a, x, b, n);*/

  /****************************************************************************
  ****************************************************************************/

  printf ("%d:%g:%g\n", LLC_NAME, llc_time, llc_resi);

  /**************************************************************************
   Libeacion de memoria */
  LiberaMatriz (a);
  LiberaVector (b);
  LiberaVector (x);

  return 0;

}
